# Carga de librerías
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.3 ✓ dplyr 1.0.0
## ✓ tidyr 1.1.0 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
library(tidymodels)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidymodels 0.1.1 ──
## ✓ broom 0.7.0 ✓ recipes 0.1.13
## ✓ dials 0.0.9 ✓ rsample 0.0.8
## ✓ infer 0.5.3 ✓ tune 0.1.1
## ✓ modeldata 0.1.0 ✓ workflows 0.2.1
## ✓ parsnip 0.1.4 ✓ yardstick 0.0.7
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidymodels_conflicts() ──
## x scales::discard() masks purrr::discard()
## x dplyr::filter() masks stats::filter()
## x recipes::fixed() masks stringr::fixed()
## x dplyr::lag() masks stats::lag()
## x yardstick::spec() masks readr::spec()
## x recipes::step() masks stats::step()
library(epiDisplay)
## Loading required package: foreign
## Loading required package: survival
## Loading required package: MASS
##
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
##
## select
## Loading required package: nnet
##
## Attaching package: 'epiDisplay'
## The following object is masked from 'package:yardstick':
##
## kap
## The following object is masked from 'package:scales':
##
## alpha
## The following object is masked from 'package:ggplot2':
##
## alpha
library(GGally)
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library(corrr)
library(knitr)
library(ggrepel)
library(kableExtra)
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
library(here)
## here() starts at /home/daro/cdatos/TPs/EEA/TP2
El dataset de precios de inmuebles proviene de Properati. El mismo ya fue filtrado por los docentes y a su vez se encuentra particionado en subconjuntos de training y test.
train_ds <- read_csv(here("/ds/ar_properties_train.csv"))
## Parsed with column specification:
## cols(
## id = col_character(),
## l3 = col_character(),
## rooms = col_double(),
## bathrooms = col_double(),
## surface_total = col_double(),
## surface_covered = col_double(),
## price = col_double(),
## property_type = col_character()
## )
head(train_ds)
dim_desc(train_ds)
## [1] "[32,132 x 8]"
summary(train_ds)
## id l3 rooms bathrooms
## Length:32132 Length:32132 Min. :1.000 Min. :1.000
## Class :character Class :character 1st Qu.:2.000 1st Qu.:1.000
## Mode :character Mode :character Median :3.000 Median :1.000
## Mean :2.722 Mean :1.427
## 3rd Qu.:3.000 3rd Qu.:2.000
## Max. :8.000 Max. :5.000
## surface_total surface_covered price property_type
## Min. : 28 Min. : 26.00 Min. : 69500 Length:32132
## 1st Qu.: 46 1st Qu.: 42.00 1st Qu.:120000 Class :character
## Median : 65 Median : 58.00 Median :169000 Mode :character
## Mean : 80 Mean : 70.61 Mean :214991
## 3rd Qu.: 99 3rd Qu.: 85.00 3rd Qu.:260000
## Max. :320 Max. :265.00 Max. :950000
unique(train_ds$l3)
## [1] "Almagro" "Palermo" "Caballito"
## [4] "Villa Crespo" "Villa Urquiza" "Barracas"
## [7] "Nuñez" "Belgrano" "Floresta"
## [10] "Recoleta" "Colegiales" "Parque Chas"
## [13] "Barrio Norte" "Villa Devoto" "Villa Ortuzar"
## [16] "Velez Sarsfield" "Parque Chacabuco" "Villa Pueyrredón"
## [19] "Villa Real" "Once" "Flores"
## [22] "San Telmo" "Las Cañitas" "Villa Santa Rita"
## [25] "Villa del Parque" "Retiro" "Parque Centenario"
## [28] "Chacarita" "Abasto" "San Cristobal"
## [31] "Liniers" "Balvanera" "Boedo"
## [34] "Villa General Mitre" "Saavedra" "Parque Patricios"
## [37] "Boca" "Paternal" "Monserrat"
## [40] "San Nicolás" "Parque Avellaneda" "Centro / Microcentro"
## [43] "Puerto Madero" "Congreso" "Monte Castro"
## [46] "Mataderos" "Coghlan" "Versalles"
## [49] "Villa Lugano" "Catalinas" "Villa Luro"
## [52] "Constitución" "Pompeya" "Tribunales"
## [55] "Agronomía" "Villa Soldati" "Villa Riachuelo"
unique(train_ds$property_type)
## [1] "Departamento" "PH" "Casa"
apply(train_ds, 2, function(x) any(is.na(x)))
## id l3 rooms bathrooms surface_total
## FALSE FALSE FALSE FALSE FALSE
## surface_covered price property_type
## FALSE FALSE FALSE
tab1(train_ds$property_type, sort.group = "decreasing", cum.percent = TRUE)
## train_ds$property_type :
## Frequency Percent Cum. percent
## Departamento 28203 87.8 87.8
## PH 3120 9.7 97.5
## Casa 809 2.5 100.0
## Total 32132 100.0 100.0
tab1(train_ds$l3, sort.group = "decreasing", cum.percent = TRUE)
## train_ds$l3 :
## Frequency Percent Cum. percent
## Palermo 4890 15.2 15.2
## Belgrano 2736 8.5 23.7
## Almagro 2486 7.7 31.5
## Caballito 2393 7.4 38.9
## Villa Crespo 2095 6.5 45.4
## Recoleta 1824 5.7 51.1
## Barrio Norte 1359 4.2 55.3
## Villa Urquiza 1319 4.1 59.4
## Flores 916 2.9 62.3
## Nuñez 887 2.8 65.1
## Balvanera 883 2.7 67.8
## Villa Devoto 527 1.6 69.4
## Colegiales 523 1.6 71.1
## Villa del Parque 500 1.6 72.6
## San Telmo 488 1.5 74.2
## Saavedra 433 1.3 75.5
## San Cristobal 417 1.3 76.8
## Parque Centenario 404 1.3 78.1
## Floresta 366 1.1 79.2
## Puerto Madero 357 1.1 80.3
## Boedo 336 1.0 81.3
## Paternal 327 1.0 82.4
## San Nicolás 325 1.0 83.4
## Monserrat 310 1.0 84.3
## Chacarita 292 0.9 85.3
## Retiro 291 0.9 86.2
## Villa Pueyrredón 285 0.9 87.0
## Parque Chacabuco 282 0.9 87.9
## Barracas 280 0.9 88.8
## Villa Luro 264 0.8 89.6
## Liniers 260 0.8 90.4
## Mataderos 259 0.8 91.2
## Once 255 0.8 92.0
## Congreso 246 0.8 92.8
## Coghlan 240 0.7 93.5
## Las Cañitas 187 0.6 94.1
## Abasto 179 0.6 94.7
## Parque Patricios 171 0.5 95.2
## Monte Castro 168 0.5 95.7
## Constitución 139 0.4 96.2
## Villa Santa Rita 128 0.4 96.6
## Villa Lugano 113 0.4 96.9
## Versalles 113 0.4 97.3
## Villa Ortuzar 109 0.3 97.6
## Centro / Microcentro 105 0.3 97.9
## Boca 104 0.3 98.3
## Villa General Mitre 99 0.3 98.6
## Parque Chas 78 0.2 98.8
## Parque Avellaneda 78 0.2 99.0
## Pompeya 65 0.2 99.2
## Velez Sarsfield 56 0.2 99.4
## Villa Real 55 0.2 99.6
## Tribunales 54 0.2 99.8
## Agronomía 51 0.2 99.9
## Villa Riachuelo 13 0.0 100.0
## Villa Soldati 9 0.0 100.0
## Catalinas 3 0.0 100.0
## Total 32132 100.0 100.0
lm_all <- lm(price ~ rooms + bathrooms + surface_total + surface_covered + property_type + l3, data = train_ds)
tidy_lm_all <- tidy(lm_all, conf.int = TRUE)
tidy_lm_all
(intercept) (-109143.7839): es una inferencia teórica del modelo, no representa la realidad, sería una propiedad sin superficie, habitaciones, barrio, etc.
rooms (-4359.7289): Al ser negativo implica una quita promedio de ese valor en el precio de la propiedad al aumentar en 1 la cantidad de ambientes.
bathrooms (34367.2123): Indica en cuanto aumenta en promedio el valor de la propiedad por agregar un baño.
surface_total (890.5585): Indica en cuanto aumenta en promedio el valor de la propiedad por aumentas en 1 m2 la superficie total.
surface_covered (1497.9310): Indica en cuanto aumenta en promedio el valor de la propiedad por aumentas en 1 m2 la superficie cubierta.
property_typeDepartamento (91485.9410): El valor de β (91485.9410) indica cuanto aumenta la función de respuesta para un departamento en comparación de una casa.
property_typePH (46220.0449): El valor de β (46220.0449) indica cuanto aumenta la función de respuesta para un PH en comparación de una casa.
l3{Barrio} : Son 56 coeficientes que indican en cuanto aumenta o se disminuye en promedio el precio de una propiedad dependiendo del barrio en donde esta ubicada.
Observación: para cada coeficiente se asume que se mantienen constantes el restos de las covariables.
Analizando los p-valores de las variables dummies asociadas a property_type estas resultan estadísticamente significativas. Por otro lado, las variables dummies asociados a la variable l3 se reparten entre algunas significativas y otras no.
tidy(anova(lm_all))
Analizando la significatividad global de las variables l3 y property_type y observando que tienen un p-valor menor a 0.5 se puede concluir que ambas son signigicativas.
glance(lm_all)
ggplot(tidy_lm_all,
aes(estimate, term, xmin = conf.low, xmax = conf.high, height = 1)) +
geom_point(color = "blue", size=3) +
geom_vline(xintercept = 0, lty = 2, color = "black") +
geom_errorbarh(color = "red", size=1) +
theme_bw() +
labs(y = "Coeficientes β", x = "Estimación")
El r.squared en el resúmen global representa el R cuadrado o coeficiente de determinación que refleja la bondad del ajuste del modelo a la variable precio. Al ser un valor cercano a 1 nos está diciendo que el modelo lm_all ajusta muy bien dicha variable.
Por otro lado, el gráfico nos permite visualizar con mayor claridad los intervalos de confiaza de de los coeficientes de las variables dummies de l3. Muchas de ellas incluyen el cero, ergo estas variables podrían considerarse no significativas.
lm_wo_l3 <- lm(price ~ rooms + bathrooms + surface_total + surface_covered + property_type, data = train_ds)
tidy_lm_wo_l3 <- tidy(lm_wo_l3, conf.int = TRUE)
tidy_lm_wo_l3
tidy(
anova(lm_wo_l3)
)
Al quitar la variable l3 la cantidad de coeficientes se reduce considerablemente, no necesitamos graficar para indentificar que ningúno de los intervalos de confianza incluye el 0. Por otro lado todos los p-values son menores a 0.5 ergo las variables se pueden considerar significativas, esto tambien se ve confirmado por el resultado del Test F.
models = list(lm_all = lm_all, lm_wo_l3 = lm_wo_l3)
purrr::map_df(models, broom::glance, .id = "model")
El R cuadrado de lm_all es superior al del modelo lm_wo_l3 por lo tanto el modelo que contiene todas las covariables es el que mejor explica la variable la variabilidad del precio.
En el análisis anterior se observó que l3 no es un buen agrupador, algunas de sus dummies eran significativas pero muchas otras no. Vamos a proceder a crear una nueva variable neighborhood que agrupe a las propiedad teniendo en cuenta el precio del metro cuadrado.
Es necesario tambien agregar otra variable auxiliar que contenta el precio por metro para luego poder hacer el agrupamiento necesario:
La valores posibles de neighborhood van a ser:
train_ds = train_ds %>% mutate(
price_m_sq = price / surface_total,
neighborhood = factor(
case_when(
price_m_sq <= quantile(price_m_sq)[2] ~ "low_price",
price_m_sq > quantile(price_m_sq)[2] & price_m_sq <= quantile(price_m_sq)[4] ~ "medium_price",
price_m_sq > quantile(price_m_sq)[4] ~ "high_price"
)
)
)
train_ds
lm_neighborhoods <- lm(price ~ rooms + bathrooms + surface_total + surface_covered + property_type + neighborhood, data = train_ds)
tify_lm_neighborhoods <- tidy(lm_neighborhoods, conf.int = TRUE)
tify_lm_neighborhoods
models = list( lm_all = lm_all, lm_wo_l3 = lm_wo_l3, lm_neighborhoods = lm_neighborhoods )
purrr::map_df(models, broom::glance, .id = "model") %>% arrange(adj.r.squared)
El modelo lm_neighborhoods supera al lm_all en el R cuadrado y por lo tanto tiene mayor explicabilidad sobre la variabilidad del precio de las propiedades.
Dada la fuerte correlación entre las variables surface_total y surface_covered, vamos a generar una nueva variable llamada sup_descubierta que representa la diferencia entre la superficie total y la cubierta.
train_ds = train_ds %>%
mutate(sup_descubierta = surface_total - surface_covered)
train_ds
Ahora con nuestra nueva variable, vamos a generar un nuevo modelo que incluya la variable sup_descubierta en reemplazo de la variable surface_total.
modelo_sup_descubierta <- lm(
price ~ rooms + bathrooms + sup_descubierta + surface_covered + property_type + price_m_sq + neighborhood ,
data = train_ds
)
tidy_msd <- tidy(modelo_sup_descubierta, conf.int = TRUE)
tidy_msd
Podemos observar como los p-valores de todos nuestros coeficientes resultan significativos en nuestro nuevo modelo.
tidy(
anova(modelo_sup_descubierta)
)
Y al ejecutar el Test F, validamos el mismo.
Ahora comparamos todos los modelos que generamos.
modelos = list(
lm_all = lm_all,
lm_wo_l3 = lm_wo_l3,
lm_neighborhoods = lm_neighborhoods,
modelo_sup_descubierta = modelo_sup_descubierta
)
purrr::map_df(modelos, broom::glance, .id = "model") %>%
arrange(adj.r.squared)
Podemos observar como nuestro nuevo modelo modelo_sup_descubierta es el que mejor explica la variabilidad.
Vamos a evaluar nuestro modelo lineal con la variable sup_descubierta
plot(modelo_sup_descubierta)
- Residuos vs valores predichos: parece existir una cierta estructura en los datos. La linealidad parece que no se respeta. Vemos la existencía de heterocedasticidad a medida que aumentan las predicciones. Se resaltan 2 observaciones. - Normal QQ plot: los extremos no se ajustan a la distribución teórica, tanto inferior izquierdo como superior derecho. - Residual vs leverage:Existen cuatro puntos con un leverage bastante alto. Y se identifican 2 puntos como posibles outliers. .
Vamos a generar un modelo para predecir el logaritmo natural de la variable price \[ log(price) = \beta_0 + \beta_1log(rooms) + \beta_2log(bathrooms) + \beta_3log(surface\_covered) + \beta_4property\_type + \beta_5barrio + \beta_6surface\_patio \]
Para esto vamos a generar nuestros nuevas variables
modelo_log = lm(
log(price) ~ log(rooms) + log(bathrooms) + log(surface_covered) + property_type + neighborhood + sup_descubierta,
data = train_ds
)
tidy_ml <- tidy(modelo_log, conf.int = TRUE)
tidy_ml
tidy(
anova(modelo_log)
)
De nuestro modelo confirmamos que todos los coeficientes resultan significativos.
Ahora comparamos nuestro nuevo modelo_log con nuestro modelo_sup_descubierta.
modelos = list(
modelo_log = modelo_log,
modelo_sup_descubierta = modelo_sup_descubierta
)
purrr::map_df(modelos, broom::glance, .id = "model")
Observamos que el mismo explica mejor la variabilidad \(Rˆ2\)
Ahora vamos a analizar el cumplimiento de los supuestos del modelo lineal
plot(modelo_log)
- Residuos vs valores predichos: parece existir una cierta estructura en los datos. La linealidad se respeta. Podemos observar un comportamiento heterocedástico aunque parece ser más uniforme que el del modelo_sup_descubierta. Se resaltan 3 observaciones (observaciones 3636,3709,7346). - Normal QQ plot: los extremos no se ajustan a la distribución teórica, tanto inferior izquierdo como superior derecho. Sin embargo parece que se comportan ligeramente mejor que el modelo_sup_descubierta (parecen estár mas suavizados los extremos). - Residual vs leverage:. se pueden apreciar 2 puntos con un leverage bastante alto y el grafico resalta 3 observaciones.
Vamos a generar un nuevo modelo utilizando como variables \(log(rooms)\), \(log(bathrooms)\), \(log(surface_total)\), property_type y neighborhood y otro similar a nuestro modelo_log reemplazando la variable neighborhood por la variable l3.
Generamos nuestro modelo modelo_log_surface_total que buscará inferir el resultado del \(log(price)\)
modelo_log_surface_total = lm(
log(price) ~ log(rooms) + log(bathrooms) + log(surface_total) + property_type + neighborhood,
data = train_ds
)
tidy_mlst <- tidy(modelo_log_surface_total, conf.int = TRUE)
tidy_mlst
En nuestro modelo podemos observar como el coeficiente de \(log(rooms)\) parece no ser significativo (su p-valor = 0.5) y su intervalo de confianza incluye el 0. El resto de nuestros coeficientes parecen ser significativos.
Debido a que neustros coeficientes para la variable \(log(rooms)\) parece que no son significativos, vamos a analizar si la interacción de la variable si lo es para nuestro modelo.
tidy(
anova(modelo_log_surface_total)
)
Y resulta que si lo es.
Analizamos el cumplimiento de los supuestos de nuestro modelo.
plot(modelo_log_surface_total)
Podemos ver como nuestro modelo se comporta de forma muy similar al modelo_log. Sin embargo hay una diferencia importante al analizar Residuos vs Leverage. Aca podemos observar como hay una especie de corte en el leverage entre los valores de 0.0007 y 0.0013 (aproximadamente).
Generamos nuestro modelo modelo_log_surface_total que buscará inferir el resultado del \(log(price)\)
modelo_log_l3 = lm(
log(price) ~ log(rooms) + log(bathrooms) + log(surface_covered) + property_type + sup_descubierta + l3 ,
data = train_ds
)
tidy_mll3 <- tidy(modelo_log_l3, conf.int = TRUE)
tidy_mll3
En este modelo, podemos observar el mismo problema que surgía al incorporar al analisis la variable l3. Sin embargo, el resto de nuestros coeficientes parecen ser estadisticamente significativos.
Debido a que neustros coeficientes para la variable l3 parece que no son significativos, vamos a analizar si la interacción de la variable si lo es para nuestro modelo.
tidy(
anova(modelo_log_l3)
)
Y resulta que si lo es.
Analizamos el cumplimiento de los supuestos de nuestro modelo.
plot(modelo_log_l3)
- Residuos vs valores predichos: parece existir una cierta estructura en los datos. La linealidad se respeta. Podemos observar un comportamiento heterocedástico (muy similar al modelo_log). Se resaltan 3 observaciones (observaciones 7346,7942, 28022). - Normal QQ plot: los extremos no se ajustan a la distribución teórica, tanto inferior izquierdo como superior derecho. Sin embargo parece haber cierto suavizado en el extremo superior derecho en comparación a otros modelos. - Residual vs leverage:. se pueden apreciar 3 puntos con un leverage bastante alto y el grafico resalta 3 observaciones (3636, 22246,26579).
Ahora vamos a elegir 2 de los modelos previamente desarrollados para evaluar y seleccionar, junto a nuestros 2 nuevos modelos, cual es el que mejor permite inferir el precio de una propiedad. Vamos a elegir los modelos que mejor expresan la variabilidad \(Rˆ2\). Estos son el modelo_log y modelo_sup_descubierta.
modelos = list(
modelo_log = modelo_log,
modelo_sup_descubierta = modelo_sup_descubierta,
modelo_log_surface_total = modelo_log_surface_total,
modelo_log_l3 = modelo_log_l3
)
purrr::map_df(modelos, broom::glance, .id = "model") %>%
arrange(adj.r.squared)
Podemos observar que el modelo que mejor explica la variabilidad es el modelo modelo_log_surface_total. Observamos el valor del \(Rˆ2\) ajustado.
Cargamos el dataset de testing y le aplicamos nuestras transformaciones.
propiedades_test <- read_csv(here("/ds/ar_properties_test.csv"))
## Parsed with column specification:
## cols(
## id = col_character(),
## l3 = col_character(),
## rooms = col_double(),
## bathrooms = col_double(),
## surface_total = col_double(),
## surface_covered = col_double(),
## price = col_double(),
## property_type = col_character()
## )
propiedades_test = propiedades_test %>%
mutate(
price_m_sq = price / surface_total,
neighborhood = factor(
case_when(
price_m_sq <= quantile(price_m_sq)[2] ~ "low_price",
price_m_sq > quantile(price_m_sq)[2] & price_m_sq <= quantile(price_m_sq)[4] ~ "medium_price",
price_m_sq > quantile(price_m_sq)[4] ~ "high_price"
)
),
sup_descubierta = surface_total - surface_covered
)
propiedades_test
Ahora vamos a evaluar nuestros modelos en training y evaluar sus predicciones en testing. Para esto es necesario hacer un análisis distinto para los modelos que utilizan \(log\) que para el modelo_sup_descubierta.
modelos_log = list(
modelo_log = modelo_log,
modelo_log_l3 = modelo_log_l3,
modelo_log_surface_total = modelo_log_surface_total
)
predicciones_training_log = map(.x = modelos_log, .f = augment)
map_dfr(
.x = predicciones_training_log,
.f = rmse,
truth = exp(`log(price)`),
estimate = exp(.fitted),
.id="modelo"
) %>% arrange(.estimate)
eval_train_sd = augment(modelo_sup_descubierta)
rmse(
data = eval_train_sd,
truth = price,
estimate = .fitted
)
Analizando los RMSE de nuestros modelos en training, podemos observar que el modelo con el menor valor es nuestro modelo_sup_descubierta seguido (en orden) por los modelos modelo_log_surface_total, modelo_log, modelo_log_l3.
Ahora vamos a evaluar nuestros modelos en testing.
# Aplicamos la función augment a los 4 modelos con el set de testing
predicciones_testing = map(
.x = modelos_log,
.f = augment,
newdata = propiedades_test
)
# Obtenemos el RMSE para los 4 modelos
map_dfr(
.x = predicciones_testing,
.f = rmse,
truth = price,
estimate = exp(.fitted),
.id="modelo"
) %>%
arrange(.estimate)
pred_log = augment(
modelo_sup_descubierta,
newdata=propiedades_test
)
pred_log
rmse(
data = pred_log,
truth = price,
estimate = .fitted
)
Volvemos a observar como el modelo con el menor valor de RMSE es nuestro modelo_sup_descubierta.
En lineas generales no podemos afirmar que nuestros modelos cumplen los supuestos del modelo lineal.
Sin embargo, podemos utilizar nuestros modelos para predecir los valores del precio de las propiedades. Ahora bien, para elegir el mejor de los modelos propuestos, vamos a fundamentar nuestra decisión en base a la evaluación de las 2 métricas que estamos observando: - \(Rˆ2\) - RMSE
Los modelos que mejor explican la variabilidad medida por \(Rˆ2\) ajustado son el modelo_log, modelo_sup_descubierta y modelo_log_surface_total. Todos estos modelos se encuentran con valores superiores a 0.91.
Sin embargo, nosotros queremos predecir nuevos datos y para ello es importante medir el RMSE para evaluar el error en la predicción. Al analizar esta metrica, observamos un salto importante entre nuestro modelo_sup_descubierta (38974.87) y sus seguidores modelo_log_surface_total (43832.07) y modelo_log (46152.97).
Por ende seleccionaría como mi mejor modelo el modelo modelo_sup_descubierta.